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Abstract 

Numerical  errors,  in  particular  bias,  in  PDF-based  particle-mesh  methods  for  turbulence 
modeling  have  been  explored.  It  is  shown  that  bias  decreases  linearly  with  the  increase  of 
the  number  of  particles,  but  increases  with  grid  refinement.  The  fluctuations  in  mean  fields 
which  are  fed  back  into  the  coefficients  of  stochastic  differential  equations  are  attributed  to 
be  the  sources  of  bias.  The  Frozen  Coefficient  approach  has  been  proposed  and  adopted  to 
pinpoint  the  sources  of  bias  in  detail.  These  results  provide  guidelines  for  improving  the 
numerical  accuracy  of  PDF  models  for  turbulent  flows. 
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INTRODUCTION 

For  complex  turbulent  reactive  flows,  probability  density  function  (PDF)  methods  have 
been  well  developed  and  offer  great  potential  [4]  [5].  In  the  application  of  these  methods,  a 
particle-mesh  method  is  used  to  solve  a  modeled  equation  for  the  evolution  of  a  PDF,  e.g., 
the  joint  velocity-frequency-composition  PDF.  This  is  a  Monte-Carlo  method  in  which  the 
PDF  is  represented  by  a  large  set  of  particles  distributed  in  the  physical  space.  This  space  is 
divided  into  a  number  of  cells  in  order  to  estimate  the  mean  fields  such  as  the  mean  velocity 
as  a  function  of  position.  Such  a  method  has  been  implemented  in  the  PDF2DV  code  [8] 
which  is  a  Fortran  code  to  calculate  the  properties  of  statistically  two-dimensional  (plane  or 
axi-symmetric)  turbulent  reactive  flows  using  the  joint  velocity-frequency-composition  PDF 
model.  In  the  present  work,  numerical  errors,  in  particular  bias,  in  PDF2DV  are  investigated. 

The  numerical  accuracy  and  convergence  of  PDF  methods  have  been  discussed  by  Pope  [6] 
and  Welton  et.al.  [9].  Usually,  weak  convergence,  i.e.,  the  convergence  of  expectations  instead 
of  the  PDF,  is  sought  for  PDF  methods.  Four  different  types  of  numerical  errors  have  been 
identified  by  considering  estimating  a  mean  quantity:  statistical  error,  spatial  discretization 
error,  temporal  discretization  error  and  bias.  The  convergence  of  numerical  solutions  to  the 
modeled  equations,  which  are  stochastic  differential  equations  (SDE),  requires  that  numerical 
errors  vanish  as  the  particle  number  per  cell  N  tends  to  infinity,  and  as  the  time  step  At  and 
grid  size  h  approach  to  zero.  It  has  been  shown  that  statistical  error  (SE)  scales  as  1/iV1/2 
[6].  It  is  reasonable  to  postulate  that  temporal  error  and  spatial  error,  behave  as  that  in 
finite  difference  method,  i.e.,  both  vanish  as  At  and  h  tend  to  zero  for  fixed  N. 

However,  in  early  experiences  with  the  PDF2DV  code  it  has  been  observed  that  there 
appeared  to  be  relatively  large  bias.  The  bias  is  the  deterministic  error  resulting  from  using 
a  finite  number  of  particles.  Its  features  in  the  PDF  methods  are  not  very  clear  yet.  This 
study  is  devoted  to  understanding  the  behavior  and  the  sources  of  bias  in  PDF2DV,  which 
will  provide  one  of  the  guidelines  for  improving  the  accuracy  of  the  code. 

The  description  of  bias  in  PDF2DV  is  presented  in  the  next  section,  and  then  the  strate¬ 
gies  for  pinpointing  the  sources  of  bias  is  described  in  the  following  section.  Two  test  cases 
are  used  to  search  for  the  sources  of  bias:  stationary  homogeneous  turbulence  and  C'ouette 
flow.  Results  for  these  cases  are  discussed.  Conclusions  are  drawn  in  the  final  section. 
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DESCRIPTION  OF  BIAS  IN  PDF2DV 

Before  presenting  the  behavior  of  bias  in  PDF2DV,  the  model  equations  and  some  numerical 
techniques  are  introduced.  A  Lagrangian  approach  is  taken  both  in  the  modeling  and  in  the 
numerical  method  which  is  a  Monte-Carlo  particle-mesh  method.  The  fluid-particle  possesses 
properties:  velocity  u  +  (t)  and  turbulence  relaxation  rate  (turbulence  frequency)  w+(t)  . 
These  properties  are  then  modelled  by  the  corresponding  stochastic  processes  u *(t)  and 
LO*(t )  based  on  the  idea  of  the  stochastic  Lagrangian  modeling  approach  [7].  (If  combustion 
is  considered,  the  stochastic  models  for  scalar  fields  are  needed.)  In  PDF2DV,  u *(t)  evolves 
according  to  the  simplified  Langevin  equation  (SLM,  [7]) 

du*(t)  =  -V(p)dt  -  Q  n  (u *(t)  -  (U))  dt  +  (Co AxQ)1/2  dW,  (1) 

where  W  is  Wiener  process,  aJ*(t )  is  solved  by  the  Ito  stochastic  differential  equation  [2] 

dco*  =  -C3  (u*  -  (u;) )  ildt  -  (w)u*SIJ,dt  +  (2C3a2nu*(io)y/2  dW.  (2) 

The  non-dimensional  source  term  S w  is  defined  as: 

S„  =  C2-  C1SijSij/(u)\  (3) 

where  .9,-,  is  the  mean  rate  of  strain.  For  all  other  terms,  model  constants  and  the  definition 
for  conditional  mean  of  turbulence  frequency  O  in  (1)  and  (2),  refer  to  [2]  [7]. 

In  PDF2DV,  the  expectation  of  a  random  variable,  e.g.  (U),  is  approximated  by  an 
ensemble  mean  estimated  by  the  cloud-in-cell  method.  Because  the  number  of  samples,  i.e., 
the  number  of  particles  per  cell  N  is  finite,  the  ensemble  mean  itself  is  a  random  variable.  Of 
course,  it  also  depends  on  time  step  At  and  grid  size  h.  Several  numerical  techniques,  e.g., 
variance  reduction  (VR)  and  time  averaging  (TAV),  are  adopted  to  reduce  the  statistical 
fluctuations  in  the  mean  fields  [8]. 

The  bias  Bq  of  a  statistics  (Q),  e.g.  (U),  is  the  deterministic  error  caused  by  N  being 
finite.  Using  (Q)  \.±u,  to  represent  the  ensemble-average  of  Q  calculated  by  finite  N,  At  and 
h.  Bq  can  be  written  as, 


B Q  (  (Q)  N,  At,/? )  (Q)  oo,At,h- 


(4) 
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Figure  1:  Time-averaged  mean  values  vs.  N  in  Couette  flow:  (a)  Mean  Velocity:  (b)  Mean 
(u>).  Grids  in  y  direction  Ng  =  21,  both  VR  and  TAV  are  off,  y/H  =  0.033. 


A  simple  analysis  by  Pope  [6]  suggests  that  bias  Bq  scales  as  N  1 
as: 


B 


so  Bq  can  be  expressed 


(5) 


where  b  may  be  a  function  of  other  factors  affecting  bias,  e.g.,  grid  size  and  other  numerical 
techniques  in  the  code. 

As  an  example,  Couette  flow  (for  flow  descriptions  refer  to  appendix  A)  is  calculated 
using  PDF2DV.  With  the  same  grid  size,  time  step  and  numerical  techniques,  the  time 
averages  (( U)n)t  and  ((lo) n)t  of  one  point,  which  are  equivalent  with  the  ensemble  mean 
of  (U)n  and  (u))i y  respectively  according  to  the  ergodic  assumption,  are  shown  in  Figure  1 
as  a  function  of  the  number  of  particles  per  cell.  Since  temporal  error  and  spatial  error  are 
independent  of  particle  number,  the  linear  relationship  in  the  figure  implies  that  bias  scales 
as  (5).  This  behavior  of  bias  is  common  in  PDF2DV.  The  scale  of  bias  as  (5)  ensures  that 
bias  approaches  zero  as  N  goes  to  infinity  and  thus  leads  to  convergence  of  the  scheme  with 
respect  to  particle  number  per  cell. 

On  the  other  hand,  the  slopes  of  the  lines  in  Figure  1  provide  the  value  of  b  in  (5)  which 
determines  the  magnitude  of  bias.  For  specific  At  and  h  (and  for  all  other  aspects  of  the 
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Figure  2:  Bias  vs.  grid  size  in  Couette  flow:  (a)  Mean  Velocity,  the  slope  is  —0.55;  (b)  Mean 
(u>),  the  slope  is  —0.80.  Both  VR  and  TAV  are  off,  y/H  =  0.033. 


numerical  method  fixed),  we  have 

(Q)nu  At,h  ~  (Q)  Iffy  At, h  =  BqNi  -  BqN2  =  b  (  —  -  —  )  '  (6) 

A  formula  for  b  is  thus  obtained 

N-jN-i 

b(At,  h)  =  — —  ((Qj^'Ath  ~  {Q)N2,At,h)  ■  (7) 

Consequently,  using  two  different  particle  numbers  per  cell,  b  and  Bq  can  be  calculated  for 
a  specific  grid  size.  By  fixing  the  time  step  and  the  numerical  techniques,  b  for  different 
grid  sizes  in  Couette  flow  is  obtained  by  (7)  and  is  plotted  against  1/Ay  in  Figure  2,  where 
Ng  is  the  number  of  cells  or  grids  in  the  domain.  Since  nonuniform  grids  are  used  in  this 
calculation,  l/Ng  is  used  here  to  represent  the  averaged  grid  size.  Figure  2  shows  that  b 
increases  with  grid  size.  This  behavior  is  observed  in  homogeneous  stationary  turbulence 
(described  in  appendix  A)  as  well.  In  this  case,  bias  is  the  deviation  of  turbulence  energy 
and  turbulence  frequency  from  the  stationary  solutions.  As  shown  in  Figure  3,  increasing 
the  number  of  cells  results  in  larger  bias. 

The  fact  that  the  bias  increases  with  grid  size  is  of  great  concern  to  PDF2DV.  If  b 
explodes  faster  with  grid  size  than  the  convergence  of  bias  due  to  increasing  N,  bias  will  not 
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Figure  3:  Bias  vs.  grid  size  in  Homogeneous  Stationary  Turbulence:  (a)  Turbulence  Energy 
A;  (b)  Mean  (uj).  Time  step  At  =  0,1s,  for  4000  steps. 

vanish  as  h  approaches  zero  and  N  goes  to  infinity.  In  the  view  of  numerical  computations, 
to  minimize  the  total  numerical  error,  one  might  want  to  increase  N(J  to  reduce  the  spatial 
discretization  error,  which  unfortunately  causes  larger  bias.  Therefore,  N  must  increase  as 
Ng  increases  in  order  to  prevent  the  bias  from  exploding  and  to  get  the  total  error  converged, 
which  thus  leads  to  unacceptably  high  computational  cost. 

STRATEGIES  FOR  EXPLORING  BIAS 

Theoretically,  the  statistical  fluctuations  in  the  coefficients  of  the  stochastic  differential  equa¬ 
tions,  e.g.,  the  Langevin  equation,  are  the  major  sources  of  bias.  The  numerical  solutions  to 
(1)  and  (2)  are  essentially  different  from  the  following  standard  problem:  given  coefficients 
a(x,t )  and  6(x,f),  an  initial  condition  A'(0)  =  xo ,  and  a  stopping  time  T  >  0,  integrate  the 
stochastic  differential  equation 

dX(t)  =  a(X(t),t)dt  +  b(X(t).l)d\V(t).  (8) 

which  has  been  well  studied  [3].  This  is  because  in  (1)  and  (2),  the  coefficients  depend  on 
the  mean  of  a  function  of  the  process,  e.g.,  (U),  which  needs  to  be  estimated  in  numerical 
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Figure  4:  Strategies  for  exploring  bias  in  PDF2DV 

computation  and  inevitably  carries  numerical  fluctuations.  The  mean  fields  in  PDF2DV  are 
calculated  by  cloud-in-cell  method  and  fed  back  into  the  coefficients.  This  feedback  causes 
bias.  The  two  questions  to  be  addressed  are: 

1.  Which  coefficients  associated  with  mean  fields  yield  bias  ? 

2.  Why  does  bias  increase  with  finer  grids  ? 

For  the  PDF2DV  code|  all  mean  fields,  terms  in  the  SDE’s  and  numerical  techniques 
which  may  be  related  to  the  above  two  issues  are  sketched  in  Figure  4. 

If,  instead,  the  coefficients  were  non-random,  independent  of  particle  properties,  the 
computed  particle  properties  at  any  time  would  be  independent,  identically  distributed, 
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and  independent  of  N.  It  follows  then  there  would  be  no  bias.  Therefore,  if  the  estimate 
of  mean  field  is  replaced  by  a  non-random  input  (say  fixed  or  frozen  (£/)),  there  is  no 
corresponding  statistical  fluctuation  in  that  mean  field  and  related  terms  in  the  SDK's  so  that 
the  contribution  of  statistical  errors  in  the  mean  field  to  bias  will  be  prevented.  According 
to  this  Frozen  Coefficient  approach,  a  mean  field  or  a  term  in  the  SDE’s  is  justified  as  a 
source  of  bias  if  freezing  it  results  in  no  bias.  The  things  to  be  tested,  which  are  possibly 
the  sources  of  bias,  and  associated  methods,  are  also  shown  in  Figure  4. 

NUMERICAL  TESTS  AND  RESULTS 

The  approach  of  Frozen  Coefficient  is  accomplished  step  by  step  through  numerical  tests. 
Because  the  fluctuations  in  the  mean  fields  are  suspected  to  be  the  major  sources  of  bias, 
first  of  all  the  calculations  are  made  by  fixing  mean  fields  in  all  coefficients  of  the  SDE’s  to 
check  the  behavior  of  bias.  If  it  disappears,  then  we  can  conclude  that  bias  arises  entirely 
from  the  fluctuations  in  the  coefficients  and  no  more  calculations  are  needed  to  search  for 
other  sources  of  bias  due  to  the  numerical  techniques.  Several  cases  are  calculated  and 
compared  to  the  base  case;  namely,  the  result  from  the  original  code  without  VR  and  TAV. 
The  conditions  for  different  cases  are  presented  in  Table  I. 

1  Sources  of  Bias:  General  Views 

In  the  calculation  of  case  1,  it  is  not  the  estimates  from  cloud-in-cell  method  that  are  used  but 
constant  values  for  mean  fields.  No  bias  is  expected  if  bias  is  only  related  to  the  statistical 
fluctuations  in  the  mean  fields.  The  result  compared  to  the  base  case  (case  0)  is  shown 
in  Table  I.  Obviously,  bias  of  velocity  and  (lo)  are  much  smaller  than  case  0  and  almost 
zero  after  fixing  mean  fields  globally.  Therefore  bias  of  these  two  variables  is  because  of 
the  fluctuations  in  the  estimates  of  mean  fields.  However,  because  the  ensemble  mean  of 
velocity  is  used  to  estimate  second  moments,  bias  of  turbulent  energy  is  increased  instead, 
which  implies  that  the  fluctuations  in  the  first  moments  are  the  source  of  bias  for  the  second 
moments.  This  will  be  discussed  further  in  other  cases. 
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2  Source  of  Bias  for  Velocity 

Fixing  Mean  Fields  for  Coefficients  in  Eq.(l) 

To  clarify  further  that  bias  of  velocity  comes  from  Eq.(l),  case  2  is  calculated  by  fixing 
all  mean  fields  in  this  equation.  As  in  Table  I,  this  case  shows  that  the  bias  of  the  mean 
velocity  is  in  the  same  order  as  case  1,  which  indicates  that  the  bias  of  velocity  is  due  to 
the  fluctuations  in  the  mean  fields  fed  back  into  the  coefficients  of  Eq.(l).  This  is  consistent 
with  the  analysis  in  the  previous  section.  To  distinguish  the  influence  on  bias  of  drift  term 
from  diffusion  term  in  Eq.(l),  the  following  cases  are  designed. 

Fixing  Mean  Velocity  for  Drift  Term 

The  mean  velocity  is  a  critical  variable  in  Eq.(l)  because  the  particle  velocity  is  forced  to 
relax  to  it,  and  the  particle  velocity,  in  turn,  is  used  to  estimate  the  mean  velocity.  This 
means  that  there  must  be  a  very  strong  interaction  between  them.  As  shown  in  Table  I 
(case  3),  the  bias  of  velocity  resulting  from  fixing  the  mean  velocity  for  the  drift  term  in 
Eq.(l)  is  very  small.  Hence,  the  fluctuation  in  mean  velocity  is  a  major  source  for  bias  of 
velocity.  One  more  interesting  observation  is  that  the  bias  of  (uj)  becomes  very  small  in  this 
case  as  well.  It  seems  that  the  behavior  of  the  bias  of  (uj)  is  dominated  by  the  mean  velocity 
although  Eq.(2)  is  the  same  in  the  form  as  Eq.(l).  This  will  be  further  discussed  later. 

It  may  be  noticed  that  in  case  3  there  is  a  huge  bias  of  k.  However,  this  phenomenon 
should  not  be  paid  too  much  attention.  When  the  drift  term  is  frozen  while  the  diffusion 
term  is  not,  in  the  equation  for  k  derived  from  Eq.(l),  — |C'oO(u*  —  (U ))dt  in  the  drift  term 
will  not  still  balance  with  the  diffusion  term,  which  may  cause  a  large  bias  of  k. 

Fixing  and  k 

In  order  to  determine  the  effect  of  fluctuations  in  the  diffusion  coefficients,  a  calculation  was 
attempted  in  which  Q  and  k  were  fixed  only  in  the  diffusion  term  in  Eq.(l).  It  turns  out  that 
the  solution  is  not  stable.  The  explanation  could  be  that  different  values  of  Q  in  the  drift 
and  diffusion  terms  give  rise  to  an  instability  of  numerical  solutions  to  the  equation.  This 
case  is  then  modified  to  let  O  take  fixed  values  both  for  the  drift  term  and  the  diffusion  term 
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and  k  be  fixed  for  the  diffusion  term  so  that  a  stationary  solution  is  obtained.  Thus,  in  this 
case  4,  k  and  Q  (but  not  (U))  are  fixed  throughout  Eq.(l).  The  bias  of  velocity  is  almost 
the  same  as  in  case  0.  The  bias  of  (lo)  is  not  much  reduced  either.  This  case  demonstrates 
that  the  diffusion  term  of  Eq.(l)  is  not  the  source  of  significant  bias  for  the  velocity  or  for 

M- 

3  Source  of  Bias  for  (a;) 

Fixing  Mean  Fields  for  Coefficients  in  Eq.(2) 

In  case  5,  the  coefficients  in  Eq.(2)  are  fixed.  The  bias  of  (u>)  is  very  small  (Table  I). 
Therefore,  the  bias  of  (lo)  stems  from  the  fluctuations  in  the  coefficients  of  Eq.(2). 

Fixing  Mean  Velocity  Gradient  for  Turbulent  Frequency  Equation 

An  analysis  given  in  Appendix  B  shows  that  bias  in  the  source  term  S w  in  the  turbulent 
frequency  model  increases  with  grid  refinement  because  of  the  fluctuations  in  the  mean 
velocity.  This  may  be  the  reason  for  the  dependence  of  bias  on  grid  size,  which  can  be 
demonstrated  by  fixing  the  mean  velocity  gradient  to  calculate  SM.  The  results  in  this 
case  (case  7)  are  compared  with  the  results  of  original  model.  Figure  5  shows  that  the 
dependence  of  bias  on  grid  size  almost  disappears  when  a  fixed  mean  velocity  gradient  is 
used  to  calculate  Sm,.  Therefore,  the  reason  that  bias  increases  when  grid  size  are  decreased 
is  that  the  fluctuations  in  the  mean  velocity  bring  about  an  additional  source  into  the 
turbulent  frequency  model  because  of  the  square  of  velocity  gradient  in  Sw.  In  addition 
to  this  observation,  it  is  shown  in  Figure  5  and  Table  I  that  bias  of  (lo)  is  almost  zero. 
Consequently,  the  bias  of  (lo)  is  dominated  by  the  fluctuations  in  the  mean  velocity. 

4  Source  of  Bias  for  Turbulent  Energy  k 

Homogeneous  Flow:  Fixing  Mean  Velocity  for  Estimate  of  Second  Moments 

In  this  case  (case  H-l  in  Table  II),  the  fixed  no-random  first  moment  (the  mean  velocity) 
is  used  to  estimate  second  moments.  Table  II  shows  that  bias  reduces  by  about  half  in 
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Figure  5:  Comparison  of  the  dependence  of  Bias  on  grid  size  in  Couette  flow:  (a)  Mean 
Velocity:  (b)  Mean  (oj);  solid  line,  base  case  (case  0);  dashed  line,  calculation  by  fixing 
velocity  gradient  (case  7).  Both  VR  and  TAV  are  off,  y /H  =  0.033. 

comparison  with  case  H-0.  As  pointed  out  before,  the  fluctuations  carried  by  first  moments 
are  the  sources  of  bias  for  A. 

Couette  Flow:  Fixing  k  for  Diffusion  Term 

This  case  is  to  test  whether  the  fluctuation  in  k  which  is  fed  back  into  the  diffusion  term  of 
Eq.(l)  is  the  source  of  bias  in  A:  or  not.  The  result  is  shown  in  Table  I  (case  6).  Apparently, 
fixing  k  for  the  diffusion  term  does  not  reduce  the  bias  of  k.  Case  4  shows,  however,  that  by 
fixing  both  Q  and  k  in  diffusion  term,  bias  of  A:  is  indeed  reduced.  Therefore,  the  diffusion 
term  is  a  source  of  bias  in  A  due  to  the  fluctuations  in  both  Q  and  A. 

Flomogeneous  Flow:  Fixing  Mean  Velocity  for  Estimate  of  Second  Moments  and 
fi,  A  for  Eq.  (1) 

The  case  H-2  is  calculated  to  clarify  further  that  the  total  bias  of  A  is  contributed  by  the 
above  sources:  the  fluctuations  in  the  first  order  moments  which  are  used  to  estimate  the 
second  moments;  the  fluctuations  in  Q  and  A  which  are  feeded  back  into  Eq.(l)  (Table  II). 
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Table  I:  Bias  in  Couette  Flow:  Test  Cases  and  Results 


Cases 

0 

1 

2 

3 

4 

5 

6 

7 

8 

Drift  Term  in  the 

Equation  for  u* 

(U) 

X 

X 

X 

0 

X 

X 

X 

X 

Diffusion  Term  in 

the  Equation  for  u* 

0 

X 

X 

X 

X 

k 

X 

X 

X 

X 

Drift  Term  in  the 

Equation  for  w* 

Q 

X 

X 

M 

X 

X 

Diffusion  Term  in 

the  Equation  for  w* 

Q 

X 

X 

H 

X 

X 

Source  Term  S b  in 

the  Equation  for  lu* 

X 

X 

X 

M 

X 

X 

Evaluation  of  k 

(U) 

Measurement 

of  Bias 

b(U) 

8.42 

0.14 

0.23 

0.37 

7.70 

3.94 

21.53 

3.95 

5.58 

bN) 

6.10 

-0.23 

1.26 

-0.96 

4.85 

0.38 

10.03 

-0.20 

4.66 

h 

0.6 

-1.96 

-2.76 

-12.0 

0.14 

0.81 

1.03 

-0.03 

-0.39 

Note.  All  calculations  are  made  on  the  conditions:  Ng  =  21,  neither  VR  nor  TAV,  y/H  = 
0.033.  6’s  are  normalized  by  mean  fields  at  the  wall,  ‘x’  denotes  the  parameter  is  fixed. 


Table  II:  Bias  in  Stationary  Homogeneous  Turbulence 


Case 

Description 

Bias  of  A: 

Bias  of  (co) 

H-0 

No  any  coefficients  fixed 

1.103 

0.0051 

H-l 

Fix  ( U )  for  estimate  of  second  moments 

0.54 

0.0071 

H-2 

Fix  ( U )  for  estimate  of  second  moments 
and  fix  Q,  k  for  Eq.(l) 

0.004 

0.010 

Note.  All  calculations  are  made  on  the  conditions:  Ng  =  3,  neither  VR  nor  TAV,  time  step 
At  =  0.1s,  for  4000  steps. 
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Figure  6:  Comparison  of  the  dependence  of  Bias  on  grid  size  in  Couette  flow:  (a)  Mean 
Velocity;  (b)  Mean  (oj);  solid  line,  base  case  (case  0);  dot-dashed  line,  calculation  by  fixing 
all  mean  fields  for  Eq.(l)  (case  2);  dotted  line,  calculation  by  fixing  Q  for  Eq.(l)  (case  8). 

5  Dependence  of  Bias  on  Grid  Size 

In  the  case  of  fixing  velocity  gradients  for  the  turbulent  frequency  model,  it  has  been  shown 
that  the  source  term  Sw  in  this  model  introduces  the  dependence  of  bias  on  grid  refinement. 
Here  two  more  cases  are  set  up  to  confirm  this  further. 

Fixing  All  Mean  Fields  for  Eq.  (1) 

It  has  been  shown  in  Table  I  that  the  bias  of  velocity  almost  disappears  in  this  case  (case  2). 
Here,  it  can  be  seen  further  from  Figure  6  that  the  bias  of  velocity  is  apparently  independent 
of  the  grid  size  while  the  bias  of  (u>)  still  increases  as  the  grid  is  refined. 

Fixing  O  for  Eq.  (1) 

By  fixing  in  Eq.(l),  the  contribution  of  the  turbulent  frequency  model  to  the  dependence 
of  bias  on  grid  refinement  can  be  further  clarified.  In  this  case  (case  8),  the  bias  of  (u>) 
still  increases  with  the  grid  refinement  (Figure  6).  In  contrast,  the  bias  of  velocity  does  not 
exhibit  such  a  behavior  anymore. 
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From  the  above  two  cases  and  the  case  of  fixing  the  velocity  gradients  for  the  turbulent 
frequency  model,  it  can  be  concluded  that  the  source  term  Sj,  causes  the  bias  of  (co)  to 
depend  on  the  grid  size  and  increase  when  the  grid  size  is  refined.  Then,  the  bias  of  velocity 
is  affected  by  (lo)  or  Q  through  Eq.(l)  so  that  it  increases  with  the  finer  grids  as  well. 

6  Summarization 

The  sources  of  bias  in  the  PDF2DV  code  have  been  summarized  in  Figure  (7)  which  provides 
a  general  picture  of  the  sources  of  bias. 

CONCLUSIONS 

In  PDF  particle-mesh  methods  for  turbulence  modeling,  three  types  of  numerical  errors  are 
identified:  statistical  error,  truncation  error,  and  bias.  The  behavior  and  the  sources  of 
bias  in  the  PDF2DV  code  applying  such  a  method  have  been  studied  in  detail  by  numerical 
experiments. 

It  has  been  verified  that  bias  is  linearly  proportional  to  Ar_1  so  that  bias  decreases  with 
increasing  particle  number,  which  is  consistent  with  analysis.  Another  observation  is  that 
bias  increases  when  the  grid  size  is  decreased.  This  is  significant  because  it  could  lead  to 
unacceptably  high  computational  cost. 

The  Frozen  Coefficient  approach  has  been  proposed  to  pinpoint  the  sources  of  bias  in 
PDF2DV.  According  to  this  approach,  the  sources  of  bias  are  found  through  fixing  or  freezing 
the  mean  fields  that  appear  in  the  stochastic  differential  equations,  i.e.  using  non-random 
value  instead  of  estimates  by  cloud-in-cell  method.  This  procedure  is  implemented  by  isolat¬ 
ing  each  term  in  the  stochastic  differential  equations.  The  source  of  bias  in  velocity  is  found 
to  be  associated  with  the  fluctuations  in  the  estimated  mean  of  velocity.  In  other  words, 
the  drift  term  in  the  Langevin  equation  for  velocity  is  the  major  source  of  bias  for  velocity. 
As  for  the  turbulent  frequency  (oj),  the  bias  is  mostly  dominated  by  the  mean  velocity  and 
its  gradient.  On  the  other  hand,  the  bias  of  second  moments  (or  turbulence  energy  k )  has 
two  sources:  the  diffusion  term  of  the  Langevin  equations  for  velocity  which  is  related  to  the 
fluctuations  in  (uj)  and  A:,  and  the  estimation  of  second  moments  based  on  the  first  moments 
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Figure  7:  The  chart  describing  the  sources  of  bias 
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which  of  course  carry  fluctuations.  Lastly,  the  dependence  of  bias  on  grid  size  has  been 
attributed  to  the  source  term  S 'w  in  the  turbulent  frequency  model.  The  analysis  also  shows 
that  the  bias  of  Sw  increases  with  grid  refinement.  This  leads  to  the  dependence  of  bias  of 
(lo)  and  thus  Q  on  grid  size,  which  furthermore  affect  the  mean  velocity  through  Eq.(l)  so 
that  it  also  depends  on  the  grid  size.  The  fact  that  bias  increases  as  the  grid  size  is  decreased 
implies  the  PDF2DV  code  is  not  unconditionally  convergent  and  is  not  acceptable. 

Discovering  the  sources  of  bias  provides  a  guideline  for  reducing  bias  and  removing  the 
dependence  of  bias  on  grid  size.  Partially  time  averaging,  variance  reduction  techniques  and 
modification  of  turbulent  frequency  model  are  possible  approaches  to  improve  the  accuracy 
of  this  PDF  particle  method.  This  will  be  discussed  in  another  report. 
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APPENDIX  A:  TEST  FLOWS 


The  test  flows  used  in  this  study  are  Couette  flow  and  stationary  homogeneous  flow.  Both 
of  them  have  the  following  features 

(i)  they  are  0-D  or  1-D  problem, 

(ii)  they  exhibit  easily  understandable  physics,  and 

(iii)  they  are  statistically  stationary. 

Because  of  these  features,  many  calculations  are  feasible  during  a  reasonable  time;  the 
sources  of  bias  are  distinguishable  from  other  issues  in  the  code  and  some  numerical  tech¬ 
niques,  e.g.  time-averaging,  can  be  tested.  The  two  flows  are  described  in  this  appendix. 

(1)  Couette  Flow 

Couette  flow  is  defined  as  the  flow  between  two  flat  plates  which  move  in  the  opposite 
direction  at  velocity  Uw  (Figure  8)  .  The  flow  is  one  dimensional.  No  pressure  gradient  in  x 
direction  is  needed  to  drive  the  flow. 

For  such  a  flow,  boundary  conditions  need  to  be  defined  on  the  wall.  Since  only  the  lower 
half  of  the  domain  is  calculated,  the  boundary  conditions  are  defined  at  the  center  line  and  at 
the  lower  wall.  In  the  frame  of  a  particle  method,  the  physical  condition  is  imposed  through 
specifying  the  particle  properties.  The  following  boundary  conditions  for  the  particle  may 
not  give  the  exact  solution  to  Couette  flow.  However,  because  we  are  mostly  interested  in 
the  numerical  features  of  PDF2DV,  we  can  put  aside  the  physical  consistency  of  calculation 
with  real  world.  The  issues  we  are  concerned  with  should  be  whether  these  conditions  give 
a  stable  and  stationary  solution. 

(a)  Center  Line 

At  the  center  line,  the  physical  conditions  are 


(U)  =  0. 


d(uv ) 
dy 


o, 


(9) 


and 

d(u) 

dy 

so  that  the  particle  conditions  are  imposed  as 


0. 


(10) 


U 


* 

R 


—  it 


* 

I  * 


(11) 
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◄ -  UL 


vR  =  (12) 

-4  =  (13) 

where  IF  denotes  particles  which  are  incident  to  the  boundary  and  ‘R’  denotes  particles 
which  are  “reflected”  from  the  boundary.  The  mean  conditions  are  obtained  as  the  average 
of  the  condition  before  and  after  reflection.  It  can  be  shown  readily  that  particle  conditions 
are  consistent  with  physical  conditions. 

(b)  Wall  Conditions 

The  boundary  conditions  at  the  wall  are  more  difficult  to  define  for  particles  than  at  the 
center  line.  A  modified  wall  function  proposed  by  Dreeben  [1]  is  used  in  this  study. 

The  conditions  for  the  mean  velocities  are 

(U)  =  Uw,  (V)  =  0.  (14) 

The  velocities  of  particles  hitting  the  wall  satisfy  the  following  conditions: 

u*R  =  iCj  +  av],  (15) 

Vr  =  -v‘n  (16) 


where  a  can  be  calculated  from  the  wall  function.  To  be  consistent  with  the  mean  conditions 
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at  the  wall,  a  must  satisfy 


—  2  (uv)t 


(v2)w  ’  (  } 

where  ‘w’  stands  for  values  at  the  wall.  Using  the  wall  function,  a  formula  for  a  can  then  be 
deduced, 

c„:{r)w/k'J- 


a  = 


kUw/u  —  In  (EXu/u) 
where  the  friction  velocity  u  is  calculated  by 


(18) 


5  =  Cl^k1'2. 

r- 


(19) 


All  other  model  constants  are  chosen  as 


C „  =  0.0841,  k  =  0.4, 
Cw  =  0.76,  E  =  8.0, 
A  =  0.05. 


(20) 

(21) 

(22) 


As  for  conditions  for  turbulence  frequency  at  the  wall,  a  modified  form  of  the  model 
proposed  by  Dreeben  et.al.  [1]  is  adopted.  The  particle  condition  is 

lor  =  (23) 

Dreeben  et.al.  deduce  the  following  formula  for  /3  [1] 

/J=-2^>+0(/A  (M) 

Here  a  simplified  form  for  /4,  i.e.  a  constant,  is  used 


/ 3  =  -0.5. 


(25) 


The  above  boundary  conditions  indeed  yield  a  stable  and  stationary  solution  to  Couette 
flow.  The  profiles  of  the  mean  velocity  and  the  mean  turbulence  frequency  from  the  above 
boundary  conditions  are  shown  in  Figure  9. 

(2)  Stationary  Homogeneous  Turbulence 

The  coefficient  +  f^o)  in  SLM  (Eq.  (1))  causes  the  turbulence  energy  to  be  dissipated 
at  the  rate  (e)  or  k(u).  If  the  |  is  omitted,  the  equation  (1)  will  pertain  to  the  hypothetical 
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Figure  9:  Calculated  profiles  of  Couette  flow:  (a)  Mean  Velocity:  (b)  Mean  (u>).  Calculation 
conditions:  Half  width  of  channel  H  =  1,  Wall  velocity  Uw  =  1,  Grids  in  y  direction  Ng  =  41 

case  of  stationary  (i.e.  non-decaying)  homogeneous  (isotropic)  turbulence.  In  this  flow,  as 
turbulence  energy  does  not  decay  and  there  does  not  exist  production  either,  theoretically 
the  solutions  to  this  flow  are  constant  against  time,  i.e.  stationary  solutions  are  expected.  To 
assure  a  stationary  solution  for  (uj)  as  well,  in  Eq.(3)  constant  C-2  is  set  to  be  zero.  Because 
the  mean  rate  of  strain  Sij  is  zero  (ideally),  a  stationary  solution  to  (uj)  is  expected.  To 
see  the  effect  of  grid  size  on  bias,  we  treat  this  flow  as  1-D  flow  instead  of  0-D,  i.e.  in  one 
direction,  say  y  direction,  multiple  grids  instead  of  one  grid  are  used. 


APPENDIX  B:  DEPENDENCE  OF  BIAS  IN  Su  ON 
GRID  SIZE 


As  in  Eq.(3),  S w  is  involved  with  the  square  of  the  strain  rate.  In  Couette  flow  and  stationary 
homogeneous  turbulence  it  reduces  to 


A 


In  (nr  V _ 

2  H  dy  )  (uo)r 


(26) 


In  the  numerical  calculation,  the  differentiation  is  replaced  by  finite  difference.  According 
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Figure  10:  Cell  and  Node  Structure 


And  then,  (U)  is  approximated  by  the  corresponding  ensemble  mean  {(  }  estimated  by 
cloud-in-cell  method.  The  bias  of  (y^)  is  defined  as 


=  //{c}„-  {a},\2\  _  gup-  (U),  \ 2 


(28) 


This  can  rewritten  as 

o  /  mi + (to,  -  muu}.  \  mi + mi  - 

\  h2  /  h2 

=  A  (<{c/}2>  -  (up)  +  i  («p>;>  -  (up  - 1  ({ u)„{u)t  -  (u)„(u),) , 

=  ^  [Wr({t7}„)  +  Var({U}.)\  ~  Cov({U}n .  {£/},), 

=  p  [  ATTT)  -  ATTT)]2  +  J?  ['''ar({£7}„)Far({£/},)]1/2 (1  -  p), 

>  |yar({J7}„)Far({J7},)]1/2®-p)|  (29) 


where  p  is  the  correlation  coefficient  between  {U}n  and  {U}s.  In  most  of  the  cases,  p  is  less 
than  one,  so  bias  of  or  Sw  is  not  zero  and  determined  by  the  variance  of  the  estimate 
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of  the  mean  velocity  and  the  grid  size.  As  far  as  there  exists  fluctuations  in  the  estimate 
of  the  mean  velocity,  Bu  y  increases  with  grid  size  refinement.  Consequently,  the  bias  of 
depends  on  the  grid  size. 
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